function  [xhat,mu,gammax,xmu,ivar,sqrtvar,var,eps,Y]=batchnorm_forward(x,gamma,beta,eps)
[N, D]=size(x);
% 1. calculate mean
mu=1./N*sum(x,1);
% 2. substract mean vector
xmu=bsxfun(@minus,x,mu);
% 3. denominator
sq=xmu.^2;
% 4. calculate variance
var=1./N*sum(sq, 1);
% 5. add eps for numerical stability
sqrtvar=sqrt(var+eps);
% 6. invert sqrtvar
ivar=1./sqrtvar;
% 7. execute normalization
xhat=bsxfun(@times,xmu,ivar);
% 8. transformation steps
gammax=bsxfun(@times, gamma, xhat);
% 9.
Y=bsxfun(@plus, gammax, beta);
% 10. store intermediate
end